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ABSTRACT 

It has been suggested that the wiggle instability (WI) of spiral shocks in a galactic 
disk is responsible for the formation of gaseous feathers observed in grand-design spiral 
galaxies. We perform both a linear stability analysis and numerical simulations to inves¬ 
tigate the effect of magnetic fields on the WI. The disk is assumed to be infinitesimally- 
thin, isothermal, and non-self-gravitating. We control the strengths of magnetic fields 
and spiral-arnr forcing using the dimensionless parameters ft and T, respectively. By 
solving the perturbation equations as a boundary-eigenvalue problem, we obtain disper¬ 
sion relations of the WI for various values of f3 = l-oo and J- = 5% and 10%. We find 
that the WI arising from the accumulation of potential vorticity at disturbed shocks is 
suppressed, albeit not completely, by magnetic fields. The stabilizing effect of magnetic 
fields is not from the perturbed fields but from the unperturbed fields that reduce the 
density compression factor in the background shocks. When T = 5% and (3 < 10 or 
J- = 10% and f3 ~ 5-10, the most unstable mode has a wavelength of ~ 0.1-0.2 times 
the arm-to-arm separation, which appears consistent with a mean spacing of observed 
feathers. 

Subject headings: galaxies: ISM - galaxies: kinematics and dynamics - galaxies: spiral 
- galaxies: structure - magnetohydrodynamics - instabilities — ISM: general - shock 
waves - stars: formation 


1. Introduction 

Stellar spiral arms in disk galaxies greatly affect galaxy evolution in various ways (e.g., Buta 
& Combes 1996; Kormendy & Kennicutt 2004; Buta 2013; Sellwood 2014). They not only trigger 
and/or organize star formation by compressing gas into spiral shocks (e.g., Roberts 1969; Roberts 
& Yuan 1970; Shu et al. 1972, 1973) but also drive secular galaxy evolution (e.g., Lin &; Shu 
1964, 1966; Toomre 1964; Elmegreen 1995; Bertin & Lin 1996; Foyle et al. 2010). Observations 
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commonly indicate that grand-design spiral arms abound with secondary structures called stellar 
“spurs” or gaseous/dust “feathers” that stretch out almost perpendicularly from the main arms 
into the interarm regions in a trailing configuration (e.g., Baade 1963; Lynds 1970; Weaver 1970; 
Elmegreen & Elmegreen 1983; Scoville & Rector 2001; Scoville et al. 2001; Kennicutt 2004; Willner 
et al. 2004; La Vigne et al. 2006; Corder et al. 2008; Silva-Villa &; Larsen 2012; Schinnerer et 
al. 2013). While spurs and feathers have similar pitch angles and thus are thought to share the 
common origin (Elmegreen 1980; see also Puerari et al. 2014), what actually forms them has been 
a matter of considerable debate. 

Existing theories for the formation of secondary structures differ in the relative role played by 
gaseous self-gravity to other hydrodynamic processes. Balbus & Cowie (1985) and Balbus (1988) 
studied a linear stability of spiral shocks to axisymmetric and non-axisymmetric modes, respec¬ 
tively, and showed that these arm substructures may form as a consequence of swing amplification 
of local density perturbations. Elmegreen (1994) showed that local perturbations can grow faster 
in the presence of magnetic fields that can remove the constraint of angular momentum conser¬ 
vation. These linear-theory predictions were confirmed by numerical simulations (Kim & Ostriker 
2002, 2006; Shetty & Ostriker 2006) that showed that both magnetic fields and self-gravity are 
essential to form secondary structures resembling feathers in the nonlinear stage. Recently, Lee 
& Shu (2012) treated a stability of self-gravitating, magnetized spiral shock as a global, rather 
than local, problem in the direction perpendicular to the arm. They neglected galactic shear for 
perturbations while keeping non-inertial terms, and found semi-analytically that such shocks are 
prone to feather formation. Although these feather-forming instabilities are termed differently as 
the azimuthal instability, magneto-Jeans instability, and feathering instability by Elmegreen (1994), 
Kim &; Ostriker (2002), and Lee & Shu (2012), respectively, they refer to the same gravitational 
instability of a rotating medium in which embedded magnetic fields play a destabilizing role. More 
recently, Lee (2014) extended the work of Lee & Shu (2012) to explore the parameter space of the 
feathering instability by varying the strength of magnetic fields and self-gravity. 

In contrast, other studies argued that gaseous self-gravity may not be indispensable for the 
formation of feathers (e.g., Johns & Nelson 1986; Wada & Koda 2004; Dobbs & Bonnell 2006, 2007). 
These authors showed that secondary structures develop even in the absence of self-gravity, with 
their spacings smaller than those resulting from the self-gravitating instabilities mentioned above. 
In particular, Wada &; Koda (2004) observed that small clumps form in strongly shocked layers 
that may grow into feathers, and termed the clump-forming mechanism the wiggle instability (WI). 
The WI appears prevailing in recent numerical simulations of spiral galaxies (Kim & Kim 2014) 
and even barred galaxies with strong dust-lane shocks (Kim et al. 2012a,b; Kim & Stone 2012; 
Seo & Kim 2013, 2014). Based on the Richardson number criterion (e.g, Chandrasekhar 1961), 
they further proposed that the WI is originated from the Kelvin-Helmholtz instability (KHI) of a 
shear layer behind the shock. However, the linear analysis of Dwarkadas & Balbus (1996) showed 
that postshock flows in the presence of rapidly varying shear are stable to the KHI. Interestingly, 
Hanawa & Kikuchi (2012) argued that the WI may result spuriously from the difficulty in resolving 
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a shock inclined to numerical grids. 

To address the nature of the WI, Kim et al. (2014, hereafter Paper I) adopted the method of Lee 
& Shu (2012) and performed a linear stability analysis of non-self-gravitating, unmagnetized spiral 
shocks. Paper I found that the WI is physical rather than numerical, arising from the generation 
of potential vorticity (PV) at a distorted shock front, known as Crocco’s theorem. Since gas in 
galaxy rotation passes through spiral shocks twice in every rotation for two-armed spirals, the PV 
contained in entropy-vortex waves can increase successively, in a Lagrangian sense, in the course 
of galaxy rotation. This sets up an Eulerian overstable normal mode along the azimuthal direction 
that grows exponentially. This is quite distinct from the KHI that relies on shear in the background 
flows. By relying on the growth of incompressible entropy-vortex modes, the WI is also different 
from the nragneto-Jeans or feathering instability mentioned above that utilizes compressive acoustic 
modes. Paper I also confirmed that the growth rate and wavelength of the most unstable mode 
found by the linear stability analysis are consistent with the results of direct numerical simulations. 

While the results of Paper I are informative to understand the nature of the WI, they are based 
on the models that do not consider magnetic fields pervasive in the interstellar medium (ISM). It 
has been well recognized that disk galaxies have a large-scale, ordered component as well as a 
small-scale, random component of magnetic fields (e.g., Wielebinski & Krause 1993; Heiles 1995; 
Beck et al. 1996; Beck 2001). In most disk galaxies, magnetic field directions based on polarized 
synchrotron radiation follow optical spiral structures fairly well, with their strength and pitch angle 
ranging typically ~ 4 —20/uG and ~ 8° —37°, respectively (see e.g., Neininger 1992; Beck et al. 1996; 
Van Eck et al. 2015). Although there is an indication that dust lanes upstream of optical spiral 
arms have strongest magnetic fields, some regular fields also extend well into interarm regions. It 
appears that small scale activities such as star formation and supernova explosions inside spiral 
arms disrupt regular magnetic fields and turn them into turbulent ones. When compressed by 
spiral shocks, these turbulent fields become as strong as, sometimes even stronger than, the regular 
component inside the arms (e.g., Fletcher et al. 2011; Houde et al. 2013; Shneider et al. 2014). For 
gas surface densities of E ~ 10 — 100 M 0 pc~ 2 inside the arms, these fields are close to equipartition 
strength with the thermal and turbulence energies, indicating that they are important in the gas 
dynamics associated with the arms (e.g., Beck et al. 1996; Chyzy et al. 2000; Beck 2007). 

In this paper, we extend Paper I by including the effects of magnetic fields on the stability of 
spiral shocks. As in Paper I, we consider a local shearing sheet of an infinitesimally-thin gaseous 
disk that is assumed to be isothermal and non-self-gravitating. To simplify the situation further, 
we consider only the regular component of magnetic fields that is initially parallel to the imposed 
stellar spiral arms. One-dimensional (ID) solutions of equilibrium magnetized spiral shocks were 
already obtained by several authors (e.g., Roberts & Yuan 1970; Kim & Ostriker 2002; Lee 2014). 
By imposing two-dimensional (2D) perturbations onto such shocks, we will show that the presence 
of magnetic fields substantially suppresses the development of the WI mainly by reducing the 
density compression factor of a background spiral shock. The tension and pressure forces from the 
perturbed magnetic fields will turn out to have a minor effect on the WI. The stabilizing effect of 
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magnetic fields is qualitatively consistent with the results of the previous magnetohydrodynamic 
(MHD) simulations (Kim & Ostriker 2006; Shetty & Ostriker 2006; Dobbs & Price 2008). We 
also run local MHD simulations to confirm the results of our linear stability analysis, and compare 
the pitch angles of the nonlinear structures formed in the simulations with the values reported by 
Puerari et al. (2014) for observed feathers in the grand-design spiral galaxies M51 and M81. 

This paper is organized as follows. In Section 2, we introduce the basic MHD equations and our 
choices of the parameters, as well as the equilibrium shock solutions we use as a background state of 
the WI. In Section 3, we derive the linearized perturbation equations for our normal-mode stability 
analysis and present the shock jump conditions. We also give the computational method to obtain 
the dispersion relations. In Section 4, we present the resulting dispersion relations for ID and 2D 
modes as well as the results of the MHD simulations, and discuss the effect of magnetic fields on 
the WI. In Section 5, we summarize our findings and discuss their astronomical implications. 


2. Steady Spiral Shocks 
2.1. Basic Equations 

We consider magnetized spiral shocks in a galactic disk and study their stability to small- 
amplitude perturbations. The disk is assumed to be infinitesimally thin, isothermal with the speed 
of sound c s , and non-self-gravitating: the effect of self-gravity will be studied in a separate work. 

Following Roberts (1969), we employ a local Cartesian coordinate system (x,y) with the x- 
and y -axes representing the directions perpendicular and parallel to a local segment of a spiral arm, 
respectively (see also Shu et al. 1973; Balbus 1988). We place the local frame at the galactocentric 
distance R and let it corotate with the arm at the pattern speed D p . Let i denote the pitch angle 
of the arm segment. By taking the local approximation (|x|, |y| <C R) and assuming that the arm 
is tightly wound (sin i <C 1), the background gas velocity in the local frame arising from galaxy 
rotation is given by 

v c = u c x + v c y = i?(D — D p ) sin i x + [i?(D — D p ) — g c Dx] y, (1) 

where D is the rotational angular velocity at R and q c = —din D/d Ini? is the local shear rate in 
the absence of the arm: the term involving q c handles differential rotation in our local models. 
The corresponding epicycle frequency is k = (R _3 d[D 2 I? 4 ]/di?) 1//2 = (4 — 2</ c ) 1 / 2 D. Note that v c 
in Equation (1) depends only on x and is independent of y (see also Kim & Ostriker 2002), which 
allows to explore the behaviors of periodic waves in the y-direction. 

Since the length scales involved in gas dynamics over galactic disks are enormously large 
compared to electrical diffusion scales (Roberts & Yuan 1970; see also Shu 1992), the magnetic 
Reynolds number is much larger than unity. In this case, one can make the “frozen-in-field” 
approximation in which the ISM is assumed to be tightly coupled to the magnetic fields. Expanding 
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the compressible, ideal MHD equations in the local frame, one obtains 


— + V • (£v t ) = 0, 

<9v 9 1 

— + v T • Vv = —cfV In £ + q c £lv x y -2!lxv - V<h s + —— 
OT 47tE 

— = V x (v T x B), 


(V x B) x B, 


together with 


V • B = 0, 


(2) 

(3) 

(4) 

(5) 


where £ = f pdz is the gas surface density, v is the gas velocity induced by the arms, and vj- = v+v c 
is the total gas velocity in the local frame. In Equations (3) and (4), B = H 1 ^ 2 B 3 D, where B 3 D 
is the midplane value of the three-dimensional magnetic fields and H = E/p is the disk thickness 
that is assumed to be constant over R (e.g., Kim &; Ostriker 2001, 2002). 

In Equation (3), <b s denotes the gravitational potential of the stellar spiral arms, for which we 
take a simple sinusoidal shape: 


<h s = 4>o cos 


( 27 TX 


\ L 


( 6 ) 


constant along the y-direction. Here, <I>o(< 0) is the amplitude and L = 27ri?sinz/m is the arm- 
to-arrn distance for m-armed spirals. We consider a domain with —1/2 < x/L < 1/2, so that the 
spiral potential achieves its minimum at the domain center (x = 0). The arm strength can be 
characterized by the dimensionless parameter 


T = 


m 


( |So| \ 


(7) 


sini \R 2 Q 2 J 

which corresponds to the maximum gravitational force due to the spiral arms relative to the cen¬ 
trifugal force of the background galaxy rotation (e.g., Roberts 1969). 

In the absence of spiral arms, the gaseous disk has a uniform surface density £ c and uniform 
thickness-normalized magnetic fields B c y parallel to the arms. We quantify the strength of the 
magnetic fields using the plasma parameter 


P = 


47rc 2 E c 

B? 


4ttP c 

B 2 

n 3D,c 


( 8 ) 


where v\ c = B 2 /4:ir'E c is the Alfven speed, and P c = c 2 p c is the mean thermal pressure of the ISM 
at the galactic plane. Adopting the fiducial values P c /kB nsj 2000 - 3000 K cm ' 3 (Heiles 2001) and 
-B 3 D, c = l-4pG (Rand & Lyne 1994) in the solar neighborhood, /3 ~ 4, but we consider a range of 
/3 to study various situations with differing field strength. 


In Appendix A, we combine Equations (2)-(5) to obtain 

d \ „ B ■ V /V x B 

« +VT ' V ) 4= toE 


( 9 ) 











where 
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V x V 7 1 + 2 fi 
S 

is PV. This states that PV is, in general, not conserved in the presence of magnetic fields 
unlike in 2D unmagnetized flows where £ remains constant along a given streamline. 

Equations (2)-(8) can be completely specified by seven parameters: q c , m , sin i, 0 p /0, c s / (RO), 
J 7 , and j3. We fix to q c = 1, m = 2, sin/ = 0.1, Q p /Q = 0.5, and c s /(RO) = 0.027, and vary 
T = 5%-10% and (3 = l-oo for our presentation below. 
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2.2. Equilibrium Shock Profiles 

We now seek for ID steady solutions of Equations (2)-(5) subject to a spiral-arm forcing with 
strength J 7 , which will be used as an unperturbed equilibrium state in Section 3. Let us denote the 
steady solutions using the subscript “0” as Eo(x), v o = uo(x)x + vo(x)y, and Bo = L>o(x)y. Then, 
the steady-state conditions yield 


Eq^to = E c tt c = constant, 


( 11 ) 


and 


uto 


duo 

dx 


° 2 s dTj o , or> d$ 

- ——-b 2Duo- 

X/Q dx 


dv 0 

uto-j— = 
dx 


K 

2 D 


dx 


Uo, 


Bp dBp 

47 tEo dx 


Bpuxp = B c u c = constant. 


( 12 ) 

(13) 

(14) 


Equations (11) and (14) imply that Bq oc So, a usual relation resulting from the conservation of 
mass and magnetic flux in one dimension. Eliminating So and B 0 in favor of uto, Equation (12) 
for m = 2 reduces to 



u|o\ du T p 
utoJ dx 


2 Duq + -RD 2 T sin 



(15) 


where 


B n 


VAO = 


V^tEq 


= VA,t 


U c \ 
UTO J 


1/2 


^ • 1/2 — 1/2 
is the Alfven speed in an equilibrium configuration. Note that uao oc B 0 oc u to . 

It can be shown that the PV in the equilibrium flows is 


£0 = 


IV X Vyo + 20 


20E ’ 


(16) 


(17) 












which is constant everywhere even in magnetized spiral shocks. This results from the fact that an 
equilibrium shock satisfies V x Bo = 0 (Eq. (9)). Note that Equation (17) requires that the local 
shear rate in a steady-state spiral shock should vary as 


_ 1 dv T o , . E 0 

qo = -Q^ = 2 - {2 - q ^ 


(18) 


indicating that shear is reversed wherever Eo/S c > 2/(2 — q c ) = 2 for q c = 1 (Balbus & Cowie 
1985; Kim & Ostriker 2002). 

The equilibrium velocity profiles can be obtained by solving Equations (13) and (15) simul¬ 
taneously. We follow the method given by Shu et al. (1973) (see also Lee & Shu 2012; Paper I). 
The detailed procedure is provided in Appendix B. Figure 1 plots exemplary profiles of equilibrium 
spiral shocks for T = 5% and 10% and (3 = oo, 10, 3, and 1. A filled circle marks the magnetosonic 
point, x mp , in each shock profile. Table 1 lists the associated values of x mp , the shock position 
x s h, and the preshock and postshock surface densities Eq - and Eq + , the density compression factor 
defined by 




g+ = 
£5“ 


TO 

s+ ■ 

TO 


(19) 


as well as the preshock Mach number A4 = u S rf 0 /c s . It is apparent that magnetic fields make 
the shock weaker by providing magnetic pressure, reducing y considerably. Since Eq - and M. 
are insensitive to /3 for fixed J-, the reduction of the compression factor due to magnetic fields 
occurs primarily by making the shock front move toward the upstream direction. Note that shear 
is reversed in the regions behind the shock front where Eo/S c > 2. The degree of shear reversal 
defined as qo in Equation (18) is larger as the shock becomes stronger, which makes the structures 
that develop as a consequence of the WI more perpendicular to the arms, as will be shown in 
Section 4.3. 


3. Linear Stability Analysis 
3.1. Perturbation Equation 


We apply small-amplitude Eulerian perturbations, denoted by Ei, u±, v%, and Bi, to an 
equilibrium shock found in the preceding section. We then linearize Equations (2)— (5) to obtain 


Do 

Dt 


Ei\ dui d In E 0 dvi _ 


DoU\ duo 


Bo 


uuo 2 5 /EW 
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and 


where 
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and mi(x,y,t) is the perturbed vector potential defined through Bi = V x (miz). 


( 22 ) 

(23) 

(24) 


Since the coefficients in Equations (20)-(23) depend only on x, we may consider perturbations 
of the form 


( Si/So \ 


/ Si(x) \ 


m 

Vi 


Ui{x) 

E(x) 


exp (—iujt + ik y y), 


(25) 


\ mi/B c / \ M\ (x) 


where u and k y are the frequency and y-wavenumber of the perturbations, respectively. Equations 
(20)-(23) then reduce to 
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where 


c/Mj . « c 

uto~; — = iwdMi H-c/i, 

ax ttro 


(27) 


(28) 

(29) 


is the Doppler-shifted frequency. We take a convention that k y is a pure real number and w is a 
complex number. 
























3.2. Shock Jump Conditions 


The applied perturbations also disturb the shock front into a sinusoidal shape. Let the shape 
of the perturbed shock front be described by 


Ci (x,y,t) = Z i exp(-iut + ik y y ), (31) 

with the constant amplitude Z\. The unit vectors normal and tangential to the instantaneous shock 
front are given by h = (1, —ik y Ci) and t = (ik y Ci, 1), respectively, while the velocity of the shock 
front is v s h = (—iwCi,0) to the first order in Ci (Dwarkadas & Balbus 1996; Lee Sz Shu 2012; Paper 
I). It is then straightforward to show that the perturbations at the perturbed shock positions can 
be written as 


s(x S h + Ci) ? 

« E 0 + S 1 + Ci^, 
dx 

(32a) 

T±(x sh + Cl) - 

d'UTO ■ > 

s WTO + + Cl , +*W_dC1) 

dx 

(32b) 

T||(x s h + Cl) - 

dvTO j. 

S VT0 + V 1 + Ci dx + ikyClUTO, 

(32c) 

-S_L (^sh H” Cl) ^ 

s ikymi — ikyCiBo, 

(32d) 

B\\(x sh + (i) - 

dmi dB 0 

s Bo o + Cl , j 
ox dx 

(32e) 


where the “_L” and “||” signs denote the components perpendicular and parallel to the instantaneous 
shock front in the stationary shock frame, respectively. All the quantities in the right-hand side of 
Equation (32) are evaluated at x = x s b. 

The jump conditions at the perturbed shock location are given by 


A s (u ± E) = 0, (33a) 

(cg + uDS- ^ 11 j = 0, (33b) 

A s (Wu||-^^ = 0, (33c) 

A s (B_ lU|| - B\\v±) = 0, (33d) 

As(B ± ) = 0, (33e) 


(e.g., Shu 1992). Substituting Equation (32) into Equation (33), one can see that the zeroth-order 
terms results in Equation (B8). The first-order terms are grouped to yield 


SoutoAs (Si) + As (XWi) + iZ\uf)A s (E 0 ) = 0, 


( u T0 + c2 s \ 

V 2uto J 


As (Si) + 
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2 u^ 0 


TO- 
duTO 
dx 


= 0, 


(34) 

(35) 
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A s (Vi) - Zi 



and 


A s (Mi) — ZiA s (Bo) /B c — 0, 


(37) 


where ui s D = ojd(x s \^)- Note that Equations (33d) and (33e) give the same result, Equation (37). 

3.3. Expansion near the Magnetosonic Point 

Since the left-hand sides of Equations (26) and (27) become identically zero at the magnetosonic 
point, the right-hand sides should also vanish at x = x mp in order for regular solutions to exist. 
Let us expand the perturbation variables near x = x mp as 


Si = a 0 + ai7] + 0(r] 2 ), 

Ui = bo + birj + 0(r] 2 ), 

Vi = c 0 + ci?y + 0(r] 2 ), 

Mi = d 0 + dir) + 0(rj 2 ), 


(38a) 

(38b) 

(38c) 

(38d) 


with coefficients oo,i, &o,i 5 c o,i, and do,i- Substituting Equations (B2) and (38) into Equations 
(26)-(29) and taking zeroth- and first-order terms in rj, one can obtain a system of five linear 
equations for the coefficients. This indicates that the solutions near x = x mp can be completely 
specified by three constants ao, bo and do- Since the resulting equations are not illuminating, we 
do not present them here. 

As in Paper I, we solve Equations (26)-(29) as a boundary value problem with eigenvector 
(Si,Ui,Vi, Mi, Zi) and eigenvalue u. Since this is a linear problem, we are allowed to take the 
amplitude of one variable arbitrarily, for which we fix Re(ao) = Im(ao) = 1 at the magnetosonic 
point. We choose three trial complex values for oj, bo, and do- This specifies the values of Si, Ui, Vi, 
and Mi as well as their derivatives at x = x mp from Equation (38). We then integrate Equations 
(26)-(29) from x = .T mp both in the downstream direction up to x = x s h + L and in the upstream 
direction to x = x s h, and apply the periodic conditions to the perturbation variables. At the shock 
front, the fourth boundary condition (Eq. (37)) determines Z\, which is in turn used to check 
the other three boundary conditions. If Equations (34)-(36) are not satisfied within tolerance, we 
return to the first step and repeat the calculations by changing bo, do, and ui, one by one, until all 
the perturbed jump conditions are met. 
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4. Results 

4.1. One-dimensional Modes with k y = 0 

We first apply the method described above to ID perturbations with k y = 0. Tables 2 and 3 
list the ten lowest eigenfrequencies with differing f3 for T = 5% and T = 10%, respectively. The 
modes are numbered in the ascending order of Re(w). Strongly magnetized shocks with small (3 
possess a small number of modes. In most cases, spiral shocks have a pure decaying mode (with 
Re(w) = 0 and Im(w) < 0), a single overstable mode (with Re(w) % 0 and Im(u;) > 0), and multiple 
underdamping modes (with Im(w) < 0). The behavior of the growth rate of the overstable mode 
with f3 is not simple. In models with T = 5%, magnetic fields tend to stabilize the ID overstable 
mode, suppressing it completely when f3 < 3. In models with T = 10%, on the other hand, Im(w) 
of the overstable mode is almost independent of (3. Figure 2 compares the eigenfunction Si between 
j3 = oo and 10 cases for five odd-n modes, showing that magnetic fields do not much affect the 
shapes of the eigenfunctions. The number of nodes in S± is 2(n — 3) for n > 5 regardless of (3. 

To verify the growth rates of the overstable modes found by our stability analysis, we run 
direct MHD simulations by utilizing the Athena code (Stone et al. 2008; Stone & Gardiner 2009). 
Among various schemes implemented in Athena, we use the constrained corner transport method 
for directionally unsplit integration, the HLLE Riemann solver for flux computation (Harten et 
al. 1983; Einfeldt et al. 1991), and the piecewise linear method for spatial reconstruction. The 
simulation domain is a ID box with length L resolved by 2048 zones. We start from an initially 
uniform surface density E c and uniform magnetic fields with (3 = 1 in the local frame described in 
Section 2. We slowly turn on the spiral arm potential amplitude and make it attain full strength 
T = 10% at t = 40/D. 

Figure 3(a) plots the temporal evolution of the gas surface density at x/L = 0. Note that E 
increases with time as the arm strength increases and saturates at about t = 50/D, after which E 
exhibits a driven-oscillator behavior. The red solid lines enveloping the density fluctuations after 
saturation correspond to the growth rate of 1.1 x 10 -2 D, very close to the value given in Table 3. 
The inset clearly shows the oscillations of E over the time interval 220 < tD < 260. Figure 3(b) 
plots the Fourier-transformed power spectrum over tD = 200-400. It is peaked at some specific 
frequencies, indicated by red arrows, equal to the real parts of the eigenfrequencies listed in Table 3. 
Note the dominance of the n = 3 mode with Re(w)/D = 2.88 in the power spectrum. This confirms 
that the results of our linear stability analysis are reliable at least for the ID perturbations. 

It is uncertain what causes the ID overstability of spiral shocks. It appears that the oversta¬ 
bility arises due to a complicated interplay among various involved agents such as spiral forcing, 
epicycle motions, thermal pressure, magnetic fields, non-uniform background density and shear, 
etc. When a spiral shock is displaced from its equilibrium position, it is forced to move backward, 
but with an increased amplitude. Notwithstanding its nature, the growth time of the overstability 
amounts to ~ 4Gyr (D/26 kms -1 kpc -1 ) -1 , which is much longer than that of the 2D wiggle in- 


stability presented below. This is also longer than the expected lifetime (shorter than ~ 1 Gyr) of 
spiral arms (e.g., Sellwood & Carlberg 1984; Oh et al. 2008, 2015; Speights & Westpfahl 2011, 2012). 
Therefore, spiral shocks can be considered stable to ID perturbations for all practical purposes. 


4.2. Two-dimensional Modes k y / 0 

4-2.1. Dispersion Relations 

We now consider 2D perturbations with k y ^ 0 to explore the WI of magnetized spiral shocks. 
Figure 4 plots the dispersion relations of ten lowest-frequency eigenmodes over 0 < k y L < 60 for a 
spiral shock with T = 5% and (3 = 10. Similarly to the ID modes, these 2D modes are numbered in 
the ascending order of Re(w) at k y = 0. The solid and dashed lines correspond to Im(cj) and Re(w), 
respectively. Note that Re(w) varies almost linearly with k y , which indicates that the eigenmodes 
can be expressed as a linear superposition of entropy-vortex waves and MHD waves (Paper I). 
Note also that each mode becomes unstable (i.e., Im(o;) > 0) in a few ranges of k y , although the 
corresponding growth rates for all n are less than 0.512. This is in contrast to the unmagnetized 
cases in which Im(u;) of the n = 7 mode keeps increasing with k y (see Fig. 4 of Paper I). 

To evaluate the quantitative effects of the magnetic fields on the WI, we compare the dispersion 
relations of overstable modes with differing (3 for the T = 5% and T = 10% cases in Figures 5 and 
6, respectively. The growth rates and wavelengths of the most unstable modes depend on T as well 
as f3 quite sensitively. The WI of unmagnetized arms is dominated by a single (n = 7 for T = 5% 
and n = 4 for 7F = 10%) mode, and this holds as long as magnetic fields are relatively weak with 
f3 > 100 for T = 5% and /3 > 5 for T = 10%. Clearly, magnetic fields reduce both the maximum 
growth rate Irn(w) max and the corresponding wavenumber k v ma . v , such that Irn(w) max /12 = 1.36 
and 1.13 occurring at ky }raax L = 100.6 and 78.1 for (3 = oo and 100, respectively, when T = 5%. 
For T = 10%, these values are increased to Im(cu) max /fl = 4.71, 2.53, 1.46, and 0.90 occurring at 
ky mav L = 198.3, 90.5, 47.1, and 31.9 for (3 = oo, 100, 10, and 5, respectively. For more strongly 
magnetized arms with (3 < 10 for T = 5% and (3 < 3 for T = 10%, on the other hand, there is no 
single dominant mode, but spiral shocks are unstable to several different modes with similar growth 
rates that become smaller with decreasing (3 and T. In this case, the range of the most unstable 
wavenumbers is k yjUiax L ~ 20-50 for T = 5% and k y ma ^L ~ 5-30 for 7F = 10%, largely independent 
of (3. For arms with /3 = 1, Im(w) max /D < 0.25 and 0.20 for T = 5% and 10%, respectively. The 
reduction of the growth rates is larger at larger k y , suggesting that magnetic fields suppress the 
growth of the WI, especially for small-scale modes. 

Figure 7 plots the eigenfunctions Si, U\. V\. M±, and Hi of the most unstable modes with 
(jj/D = 84.4 + 1.13* and k y L = 78.1 for (3 = 100 in the left panels, and with cu/fl = 59.1 + 0.50* 
and k y L = 53.4 for j3 = 10 in the right panels. Also plotted as the black solid lines in the bottom 
panels are the solutions of Equation (A9), which are in good agreement with the direct numerical 
results. The spiral forcing is fixed to F = 5% for both cases. When (3 = 100, the amplitudes 
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of the eigenfunctions decrease monotonically with the distance downstream from the shock front, 
similarly to the unmagnetized cases (Paper I). For more strongly magnetized spiral shocks, however, 
the eigenfunctions do not decay monotonically. Although | Hi | decreases with x right after the shock 
front, the magnetic tension and pressure forces from the perturbed fields cause the fluctuations in 
the amplitude of Hi further downstream. 

Figure 8 overlays the configuration of the perturbed magnetic fields shown as black solid lines 
over the real perturbed PV constructed as 

Re(£i) = Re(Hi) cos (k y y) - Im(Hi) sin (k y y), (39) 

in the regions with —0.2 < x/L < 0.3 and 0 < y/\ y < 1, with \ y = ‘2'n fk y ^ for the models shown 
in Figure 7. 1 White dots in both panels trace the wavefronts of the perturbed PV, clearly showing 
a discontinuity in k x across the shock located at x/L = —0.07 and —0.08 for (3 = 100 and 10, 
respectively. Gas motions associated with entropy-vortex modes and MHD modes not only bend 
magnetic fields but also compress them in the postshock flows. This in turn breaks the conservation 
of the PV and results in a non-monotonic behavior of |Hi|. Consequently, the net reduction of the 
perturbed PV from one shock to next is smaller for shocks with stronger magnetic fields, which 
should be matched by the jump of the perturbed PV, (Hi), across the shock front. Note that 
the perturbed fields are strong inside the regions of positive Re(£i) created by the WI, and reverse 
the directions in the regions between them. 


4-2.2. Effects of Magnetic Fields on the WI 


For unmagnetized shocks, Paper I derived an analytic expression for the change of the per¬ 
turbed PV across a disturbed shock, demonstrating that the shock distortion is indeed a source of 
the PV (see also Hayes 1957; Kevlahan 1997). Appendix C derives a similar expression by combin¬ 
ing Equations (34)-(36) for the perturbed PV jump across a magnetized spiral shock. The derived 
Equation (Cll) recovers Equation (A8) of Paper I when /3 = oo. 

While A s (Hi) formally depends on all of the five perturbation variables, we find empirically 
that the U \, Vj, and Z\ terms dominate in Equation (Cll). Since \duTo/dx\ -C |wfj| for overstable 
modes, one can approximately write 


^ (so, y- z\- 

JV, - ^ ' V-.ll + .'l/M J 

Lrvnj Zjh 


where the coefficients are given as 


£u, H = 


y 2 


(40) 


(41) 


1 Note that the configurations of the perturbed field lines shown in Figure 8 are only for an illustrative purpose: 
they do not well trace the total field lines when the unperturbed component is stronger than the perturbed one. 
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U T0 

. s {n - l ) 2 
£z, H = - 2 • 




£z,M = 


2iuj s D AB 
m(m +1) 


[(/i + l) 2 -/i(/i - 1)S] 


(42) 

(43) 

(44) 

(45) 


with A and B defined by Equations (C4) and (C5), respectively. The subscripts “H” and “M” 
in the coefficients stand for the hydrodynamic and magnetic contributions to the perturbed PV, 
respectively. For unmagnetized shocks, B = 0, so that £u,m = £z, M = 0. 


Note that the U\ term in Equation (40) results from the tangential variation of the velocity 
perpendicular to the shock, which plays a stabilizing role by reducing the PV by shock compression. 
The £u ,m term explains the PV change by the perturbed magnetic pressure (Eq. (33b)). That both 
£u ,h and £jj ,m are always positive implies that the perturbed magnetic pressure stabilizes the 
WI. This can be seen more explicitly by considering a special case where the PV is retained only 
in U\, with the other perturbation variables taken zero. Then, one can show that = 

(2/i - l)/n 2 - £ um , which is always less than unity. The V\ term in Equation (40) comes simply 
from the discontinuity of the radial wavenumber across the shock, independent of the magnetic 
fields, which also stabilizes the WI (Paper I). 


On the other hand, the Z\ term originates from the shock deformation along the tangential 
direction. This is a source for the PV generation, leading to the WI. The £z,m term is due to the 
perturbed magnetic stress (Eq. (33c)). Since £z,h and £z,m always have the same sign, the stress 
of the deformed magnetic fields destabilizes the WI. Therefore, the role of the perturbed magnetic 
fields differs in the U\ and Z\ terms. As plotted in Figure 9, however, the ratios £u t M/£u,M and 
£z,m/£z,h are quite small, especially for large (3 and T. The maximum contribution of the magnetic 
terms is less than 30% of the hydrodynamic terms, which occurs when T = 5% and (3=1. This 
indicates that the effects of perturbed magnetic fields themselves to the WI are not significant. 


What is then responsible for the reduction of the growth rates in the presence of magnetic fields, 
as shown in Figures 5 and 6? It is through the magnetic pressure of the unperturbed background 
fields that tend to reduce the compression factor /r (Table 1; see also Figure 3 of Lee (2014)). The 
amount of the PV production is smaller when a shock is weaker for fixed Z\ (Eqs. (44) and (45)). 
More importantly, the amplitudes of the perturbation variables obtained by integrating Equations 
(26)-(29) decay less with x for smaller p. (e.g., Fig. 7), resulting in larger values of || and |V( S- | 
in more strongly magnetized shocks. This enhances the stabilizing role of the velocity terms relative 
to the destabilizing Z\ term in Equation (40). Consequently, spiral shocks with stronger magnetic 
fields become more stable to the WI. 
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4.3. Numerical Simulation 

To check the results of our 2D linear stability analysis, we run direct MHD simulations using 
the Athena code. As a simulation domain, we set up a rectangular box with size L x 21 that is 
resolved by a 2048 X 4096 grid. 2 Initially, we construct a ID shock profile found in Section 2.2 
as a background state. We then apply small-amplitude perturbations to the background density 
that are realized by a flat-power Gaussian random field with a standard deviation of 10” 3 So. We 
take the shearing-periodic boundary conditions at the x-boundaries and the periodic boundary 
conditions at the //-boundaries (e.g, Hawley et al. 1995; Kim &; Ostriker 2002, 2006). 

Figure 10 displays snapshots of the gas surface density in logarithmic scale and the configura¬ 
tions of the magnetic fields in the regions with —0.2 < x/L < 0.35 and 0.5 <y/L< 1.2 for models 
with /? = 100, 10, 5, and 3 at tQ = 8, 16, 25, 40 when the WI saturates, respectively. The number of 
nonlinear structures along the y-direction that develop most strongly over the simulations domain 
is 21, 15, 12, and 9, corresponding to the wavenumber of k V m^ L = 66.0, 47.1, 37.7, and 28.3 for 
/3 = 100, 10, 5, and 3, respectively. Note that the magnetic fields bend around nonlinear structures 
that are more strongly magnetized than the surrounding regions. 

Figure 11 compares the temporal evolution of the maximum surface density measured at x/L = 
0 in these models. The fastest growing mode in each model has a slope of 0.48, 0.22, 0.17, and 
0.13, corresponding to the growth rate of Irn(o;) max /D = 1.11, 0.50, 0.39, and 0.30, for /3 = 100, 
10, 5, and 3, respectively. These numerically-measured wavelengths and growth rates are marked 
as star symbols in Figure 5, in good agreement with the analytic results for the n = 7, 10, 9, 
and 7 modes, respectively. Since various overstable modes with different n have similar maximum 
growth rates, which mode the system picks up should also depend on the initial power imposed 
by specific perturbations. In the model with /3 = 3, for example, the initial density perturbations 
with k y L = 25.1 corresponding to the maximum growth rate of the n = 7 mode are about an order 
of magnitude larger than those with k y L = 37.8 corresponding to the peak of the n = 9 mode, 
emerging most strongly in the nonlinear stage, despite having a slightly (~ 6.7%) smaller growth 
rate. 


Finally, we remark on the level of turbulence generated by the WI. In the simulations de¬ 
scribed above, the density-weighted velocity dispersion, a y = (fS(v — vq) 2 dxdy/ f Edxdy) 1 / 2 , in 
the direction parallel to the arms is found to be ~ 1.4,1.1, 0.8, and 0.7kms -1 at the time when the 
WI saturates for /3 = 100, 10, 5, and 3, respectively, which are interestingly equal approximately 
to lm(uj)mav/ky mav 3 Due to nonlinear interactions and mergers of clumps created by the WI, a y 


“By also running models with size L x L, we have confirmed that the numerical results are insensitive to the 

domain size along the {/-direction since fcy.maxT -C 1. 

3 The density-weighted velocity dispersion in the direction perpendicular to the arm is not entirely due to the WI 
because of the contamination by ID overstable modes that make the shock move back and forth in the ^-direction 
(Section 4.1). 
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increases further by about a factor of 2. This suggests that the WI can be a considerable, although 
not dominant, source of turbulence energy in the ISM since the energy injection occurs to the 
densest part of the gas. 


5. Summary and Discussion 

We have investigated the WI of magnetized spiral shocks in a galactic disk using both a linear 
stability analysis and direct MHD simulations. This is a straightforward extension of Paper I that 
studied the case of unmagnetized spiral shocks. We assume that the gas disk is infinitesimally thin, 
isothermal, and non-self-gravitating. We parameterize the strengths of the stellar spiral arms and 
magnetic fields using the dimensionless parameters T and /?, respectively (Eqs. (7) and (8)). As 
background states, we first obtain the steady spiral shock profiles with differing T and (3. We then 
impose small-amplitude wiggling perturbations to the equilibrium spiral shocks, and calculate their 
dispersion relations as a boundary-value problem with eigenfrequencies. Our local MHD simulations 
readily pick up the most unstable modes for a given set up parameters, with the numerical growth 
rates very close to the results of the linear stability analysis. 

The existence of the overstable modes proves that the WI is physical in origin, resulting from 
the accumulation of the perturbed PV from a distorted shock front. Our results show that magnetic 
fields suppress the WI, but not completely, at least for (3 > 1. The stabilizing role of magnetic fields 
is not from the perturbed fields but directly from the background unperturbed fields that tend to 
reduce the shock compression factor fj, by exerting magnetic pressure. When magnetic fields are 
relatively weak, the WI is dominated by a single dominant mode. When T = 5%, the most unstable 
mode has a growth rate Im(cu) max /D = 1.36 and 1.13 occurring at k y m ^ L = 100.6 and 78.1 for 
/3 = oo and 100, respectively. When T = 10%, these values are increased to Irn(o;) max /n = 4.71, 
2.53, 1.46, and 0.90 occurring at ky jiaax L = 198.3, 90.5, 47.1, and 31.9 for (3 = oo, 100, 10, and 5, 
respectively. For more strongly magnetized arms with (3 < 10 for T = 5% and (3 < 3 for T = 10%, 
on the other hand, several overstable modes have similar growth rates that become smaller with 
decreasing (3, while the wavelength range of the most unstable modes is k y , max L ~ 20-50 for T = 5% 
and knj m^ L ~ 5-30 for T = 10%, insensitive to (3. 

We have found that magnetic fields stabilize the WI at small scales by reducing the shock 
compression factor. When /3 = 10, the most unstable wavelength is decreased by a factor of about 
2 and 4 for T = 5% and 10%, respectively, compared to the unmagnetized cases. The stabilization 
of the WI by magnetic fields has already been observed in previous nonlinear simulations of spiral 
galaxies. For instance, Kim & Ostriker (2006) ran 2D local models for the formation of feathers and 
provided the magnetic field topology at the nonlinear stage resulting from the vorticity generation 
near the shock front. They further showed that vortical clumps produced merge into massive clouds 
with mass ~ 10 7 M 0 each in the presence of self-gravity. Grid-based global simulations by Shetty 
&: Ostriker (2006) found no indication of the development of the WI in a model with (3=1 and 
T = 10%. Based on our results, this is presumably not because magnetic fields completely suppress 
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the WI but because its growth time is longer than the simulation time span (two orbits) in their 
models with very weak (~ 0.1%) initial perturbations. 1 By running particle-based simulations, 
on the other hand, Dobbs & Price (2008) found non-axisymmetric structures, albeit weak, sill 
grow in models with (3 = 1. This is probably because they might have large-amplitude density 
perturbations arising from the Poisson noises in the initial particle distributions, helping the WI 
readily manifest in their simulations. 

It is interesting to compare the predicted wavelength X v m ^ = 2n/k y m ^ from our linear 
stability analysis with the observed spacings of feathers. By analyzing the Hubble archival data, La 
Vigne et al. (2006) measured the feather spacings in grand-design spiral galaxies. They reported 
that M51 and M74 have the feather spacing that increases from ~ 0.2 kpc in the inner regions 
(R r\j lkpc) to ~ lkpc in the outer regions (R ~ 10kpc). Adopting the arm pitch angles of 
i = 21.1° for M51 (Shetty et al. 2007) and i = 15.7° for M74 (Gusev <fc Efremov 2013) (see 
also Honig &; Reid 2015), these correspond to \ V m ^/L ~ 0.1-0.2, consistent with our results for 
magnetized arms with (3 < 10 for T = 5% or with f3 ~ 5-10 for T = 10%. 4 5 Elmegreen & Elmegreen 
(1983) presented the separations of H II regions and H I superclouds along arms in various spiral 
galaxies. Taking m = 2 and i = 20° arbitrarily, the observed separations listed in their Table 2 are 
distributed in the range Ay !max /L = 0.1-1, with a mean value of ~ 0.4, a factor of ~ 2-4 larger 
than the mean feather spacings mentioned above. This may indicate that H II regions and H I 
superclouds represent highly nonlinear structures created by mergers of WI-induced feathers. 

Recently, Puerari et al. (2014) introduced a new method to locally determine the distributions 
of the pitch angles of spiral arms and their substructures. Applying the method to 8 /xm Spitzer 
images of M51 and M81, they found that the mean difference between the pitch angles of the main 
spiral arms and the interarm secondary structures is A* ~ 10°-30°. We apply the same method 
to the results of our MHD simulations to calculate the pitch angles of the nonlinear structures 
resulting from the WI off the main arms. Figure 12 plots A i when the structures saturate in 
models with (3 = 3-100 and T = 5%, averaged over the immediate postshock regions from x = x s h 
to x = x s h + L/4. Note that A i is smaller in more strongly magnetized models. This is mostly 
because of shear reversal in the postshock regions, discussed in Section 2.2, which is larger for 
stronger shocks. Since shear reversal tends to decrease k x , Ai = tan ~ 1 (k y /k x ) is smaller for weaker 
shocks with smaller (3. The numerical values of A i ~ 20°-35° for (3 ~ 3-10 are close to the values 
obtained by Puerari et al. (2014) for M51 and M81. 

As discussed above, the theoretical predictions of our analysis based on 2D, non-self-gravitating 


4 Figure 6 indicates that the maximum growth rate of the WI is Im(tu) ma x — 0.25P for /3 = 1 and J- = 10%. Thus, 
the amplitude of the perturbations that start out initially at 0.1% is expected to grow via WI to 10~ 3 exp(0.25 x 47r) ~ 
0.02 after two orbits, which is too small to be evident in the models of Shetty & Ostriker (2006). 

5 In M51, the total field strength in the arm regions is ~ 20-25/xG (Fletcher et al. 2011). Taking the mean gas 
surface density of 10 2 -10 2 ' 5 M 0 pc -2 (Meidt et al. 2013) and a disk thickness of H = 200 pc, this corresponds to 
darm ~ 0.4-1.7. 
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models are similar to the observed properties of feathers, suggesting that the WI may be responsible 
for the feather formation in the presence of magnetic fields. However, it still remains to be seen 
whether the WI would grow into observed interarm features in real disk galaxies for the following 
reasons. First, Kim &; Ostriker (2006) showed that spiral shocks in vertically-stratified disk exhibit 
non-steady flapping motions in the direction perpendicular to the arms. Together with strong 
vertical shear, these shock flapping motions tend to disrupt coherent vortical structures and thus 
prevent the growth of the WI. 

Second, Lee & Shu (2012) and Lee (2014) recently conducted a parameter study of feathering 
instability of a magnetized, self-gravitating spiral shock. In particular, Lee (2014) found that the 
most unstable mode of the feathering instability is ~ 530 pc in their M51 arm, similar to the ob¬ 
served feather spacing reported by La Vigne et al. (2006). He also found that the growth rate of the 
feathering instability is typically ~ Cl, which is larger than the growth rate (~ 0.2-0.312) of the WI 
for 0 = 1 studied in this work. Without considering the effects of the background shear, however, 
they were unable to study the combined effects of the wiggle and feathering instabilities. It is 
possible that the WI grows fast and provides seeds for the onset of the feathering instability. Or, 
the WI is completely suppressed by shock flapping motions in stratified disks. To fully understand 
the feather formation, therefore, it is necessary to study the stability of spiral shocks in a verti¬ 
cally extended, self-gravitating, magnetized disk that harbor not only the wiggle and feathering 
instabilities but also the Parker instability. 

We acknowledge helpful discussions with J. Kim and D. Ryu. We are also grateful to the referee 
for a thoughtful report. This work was supported by the National Research Foundation of Korea 
(NRF) grant, No. 2008-0060544, funded by the Korea government (MSIP). The computation of 
this work was supported by the Supercomputing Center/Korea Institute of Science and Technology 
Information with supercomputing resources including technical support (KSC-2014-C3-003). 

A. Potential Vorticity 

Here, we derive an equation for the evolution of the perturbed PV. Equation (3) can be 
decomposed as 



(Al) 



(A2) 


where 



(A3) 


is the Lagrangian time-derivative, and we have used the identity v-p • Vv = vy • Vv^ + q c CluT■ 
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Differentiating Equations (Al) and (A2) with respect to y and x, respectively, and then sub¬ 
tracting the resulting equations, we obtain 


~~ (V x V 7 1 + 2D) = (V x v T + 2D) In £ + ^ / V x B 

Dt Dt 47r 


which can be simplified to 


D£ _ B • V (V x B 
Dt 47tE 


(A4) 


(A5) 


where £ is the PV defined in Equation (10). Equation (A5) states that £ is in general not conserved 
in magnetized flows, which is in contrast to the hydrodynamic case where D£/Dt = 0. 


After linearizing Equation (10), one can write the perturbed PV as 


6 = 


|V x vi 


-£ 0 ^, 

Eo 


(A 6 ) 


where £0 is the PV in the background unperturbed shock (Eq. (17)). Analogous to Equation (25), 
we define the amplitude Hi of £1 as 


= £i(z, y ,ty^yy = ^(^- ik y u !) - fo-Si. 


So \ dx 

In terms of the perturbation variables, Equation (A5) then becomes 


d 


, - ik v c 2 s r 

_ !a , D+UT 0 _j = i = _ 


d 2 

^0 + ky ) Af] — 


J_d^o d?M\ 
T, r dx 1 dx 2 


(A7) 


(A 8 ) 


The terms in the first parentheses in the right-hand side of Equation (A 8 ) are due to magnetic 
tension, while those from magnetic pressure are given in the second parentheses. 


A formal solution of Equation (A 8 ) is given by 

Hi(x) = (Ef^ + T t + T P ) • e~ rIm(w) exp (if k XjV (x)dx\ (A9) 

V Jx sh J 

where is the PV at the immediate postshock regions, r = f x dx/itTo is the Lagrangian time 

1 •''3'sh 

elapsed from the shock front, and k x (x) = (Re(cn) — VTok y )/uTo is the local x-wavenumber. In 
Equation (A9), Tp and Tp are defined by 


ik v<% 


Puto^o Jx sh \dx' 2 


In So + k 2 y ) Ml • e T (*') Im («) 


exp — 1 


i f k x (x")d. 


lx" ) dx , (A 10 ) 


and 


ik y° 2 s 


PUTO^O Jx ah V S c dx' 


1 dSo c , d 2 Mi 
Sl + 


e r(V)Im(a,) exp 


i f k x {x")dx" j dx', (All) 

dx sh J 


representing the contributions to the perturbed PV due to magnetic tension and magnetic pressure, 
respectively. 
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B. Method of obtaining ID Equilibrium Shock Profile 


The coefficient of duTo/dx in the left-hand side of Equation (15) vanishes at the magnetosonic 
point, x = x mp , where 


ut o 


uto 


c s u c 

U T0@ ) T = r 


= 0 , 


(Bl) 


which is a cubic equation in uto • A subsonic gas should accelerate to a supersonic speed through 
the magnetosonic point. To obtain smooth solutions, we expand uto and vo around x = x mp as 


uto = a 0 + air] + a2r] 2 + 0(r] 3 ), 
vo = 7o + 7i7 + 72 r / 2 + 0 (j 7 3 ), 


(B2a) 

(B2b) 


where i] = x — x mp and ao denotes the positive real root of Equation (Bl) for given u c , c s , and /3. 
The coefficients am 2 and 70 , 1,2 are to be determined by series expansions near 77 = 0. 

Substituting Equation (B2) in Equations (13) and (15) and keeping terms up to the second 
order in 77 , one obtains 
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(B3) 

(B4) 

(B5) 

(B 6 ) 

(B7) 


With these coefficients, uto an d vo in Equation (B2) vary smoothly near the magnetosonic point. 

The equilibrium spiral-shock profiles should satisfy the following jump conditions at the shock 
front, x = x s h: 


A s (uto^o) 

= 0 , 

(B 8 a) 

As ^(c^ + u^ 0 )T, 0 + 

= 0 , 

(B 8 b) 

A s (u 0 ) 

= 0 , 

(B 8 c) 

A s ( utoBq ) 

= 0 , 

(B 8 d) 


where A s (/) = f s+ — f s , with the superscripts “s+” and “s—” representing the quantities 
evaluated at the immediate behind (x = x s h + 0 ) and ahead (x = x s h — 0 ) of the shock front, 
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respectively. Note that Equations (B8a) and (B8d) are satisfied automatically from Equations (11) 
and (14). 


Equilibrium profiles of magnetized spiral shocks can be constructed as follows. For given T and 
j3, we first choose the magnetosonic point x mp arbitrarily, and integrate Equations (13) and (15) 
starting from x = x mp in both the upstream and downstream directions. We apply the periodic 


boundary conditions at x/L = ±1/2 and determine the shock position x s h by imposing Equation 
(B8c). We then check whether Equation (B8b) is also fulfilled. If not within tolerance, we return 
to the first step to change x mp until all the jump conditions are satisfied. 


C. Jump of Potential Vorticity at the Distorted Shock Front 


Here we derive an algebraic expression for the jump condition of the perturbed PV, A s (Si) = 
S^ + — at the shock front (x = x s h) in terms of the preshock variables (Sj 5- , U {~, VJ 8- , and 
M{~) and Z\. 

With the help of Equation (19), the jump conditions (Eq. (B8)) of the background steady 
shock can be combined to yield 



(Cl) 


for the postshock Mach number, and 



(C2) 


for the preshock Mach number. For later purposes, we calculate 



(C3) 


where 



(C4) 


and 


(C5) 


(C5) 


Note that A4 = n 1 ^ 2 , A = 1/(1 ± fj,), and B = 0 for isothermal hydrodynamic shocks. Table 1 
shows that A4 is insensitive to /3 for given T. Both A and B are positive definite. 


It is straightforward to show that Equation (34) becomes 


S[ + = S{- + [Ui~ - -(fjL- 


(C6) 
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while Equation (37) is simplified to 

M s + = M{~ + (fi - l)(u c /u S rf Q )Z 1 . 

On the other hand, Equation (35) utilizing Equations (Cl) and (C7) becomes 

—= ——[/a + 1 + (^ + 2 )B\ U'l 

( ' ,+1)e C-sr+^SM?- 


I 1 


u. 


iu? D 


(At + 1 ) S 


+ (/x - 1)S 


/ \ Uj 1 dUrrir-. 

Z\ + (1 + n — B) - — — Z\. 

/r ax 


(C7) 


(C8) 


Plugging Equation (C8) into Equation (C6) allows to write S B+ in terms of the preshock quantities. 
Lastly, Equation (36) together with Equation (C7) gives 


Vl s+ = vr-(»- 1) 


K 2 U r 


— ik v 


m T 0 


U T0 _|_ u c ^ 


H PM 2 J 


+ (C9) 


Combining Equations (28) and (A7), the jump of the perturbed PV can be written as 


A S (E 1 )=--^ 


iky (U{ 


TS + 




-ur 


=wr - un 


zkyUrpQ K/ Ur 
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iq c QLk 


l T o 


20 (-u T0 )^S 0 

s_r_s_r\j^ 

ys- ys- / n s ~T s ~ y 1 1 ' 

Z^o Zjq / “TO^O 


(CIO) 


y vr- 


ni S ~Y S ~ 1 

a T 0^0 


iky 


/LM 2 £, 


^ - 
' ax 


du : 


T0 Mr 


dx 


Using Equations (C6)-(C9), we eliminate S ^ + , Vj' s+ , and M^ + in Equation (CIO) and arrange 
the terms to obtain 


a-(ho _ c sr c ur c vr e Mr e z°- 

lh,y Zjq Zjq Zj q Zjq Zjq 


with the coefficients defined as 


£s = 


£u — 


r+i)AB ( k 2 \^ s - 

r v m 2 J to ’ 

(^r +2AB ( 1+ i + (t±m), 

r \ m m / 


q c klL 

tv =- 

u T0 


(Cll) 

(C12) 

(C13) 

(C14) 


c _ 2AB u T0 

M ~ ^(i + / i )^7 


V - 1 + |) ^ + [! - M 2 + M 1 + M + 2/x 2 )tf] } , (C15) 
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0* - i) 2 
Mi + M 



In the limit of vanishing magnetic fields (.Ad 2 —> fi and A —> (1 + n) 1 ), £s = £m = 0, and 
Equation (Cll) reduces to 



(C17) 


identical to Equation (A8) of Paper I. Paper I showed that the uj s d Z\ term in Equation (C17) is 
responsible for the production of the PV at a disturbed shock, while the terms involving Ul~ and 
V^~ suppress the PV due to shock compression and background shear, respectively. 
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Table 1. Properties of Equilibrium Spiral Shocks 


T 

P 

^mp/ L 

^sh/ L 

Sg-/£ c 

^ + /£ c 

h 

M 


oo 

0.007 

-0.069 

0.546 

6.31 

11.6 

3.40 


100 

0.008 

-0.070 

0.545 

6.12 

11.2 

3.41 

0.05 

10 

0.017 

-0.078 

0.538 

5.01 

9.29 

3.45 


5 

0.030 

-0.085 

0.534 

4.34 

8.12 

3.48 


3 

0.057 

-0.093 

0.528 

3.79 

7.17 

3.52 


1 

0.193 

-0.125 

0.512 

2.62 

5.12 

3.63 


oo 

0.115 

-0.005 

0.316 

10.9 

34.5 

5.87 


100 

0.118 

-0.006 

0.316 

10.4 

32.7 

5.88 

0.10 

10 

0.143 

-0.013 

0.314 

7.80 

24.8 

5.92 


5 

0.165 

-0.019 

0.312 

6.54 

20.9 

5.94 


3 

0.187 

-0.026 

0.311 

5.58 

17.9 

5.97 


1 

0.249 

-0.054 

0.310 

3.70 

11.9 

5.99 


Note. — For the arm and galaxy parameters of q c = 1, m = 2, 
sinz = 0.1, flp/fl = 0.5, and c s /(R£l) = 0.027. The spiral potential 
is minimized at x = 0. 
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Table 2. Eigenfrequencies of One-dimensional Perturbations for T = 0.05 


mode 


(3 = oo 


(3 = 100 


o 

T— 1 

II 

Re(w)/f2 

Im(w)/f2 

Re(w)/f2 

Im(w)/fi 

Re(w)/f2 

Im(w)/fi 

1 

0.000 

-3.727 x 10" 1 

0.000 

-3.746 x 10” 1 

0.000 

-3.760 x 10” 1 

2 

0.692 

-1.782 x 10 -1 

0.695 

-1.799 x 10" 1 

0.714 

-1.859 x 10” 1 

3 

1.496 

+1.222 x 10 -3 

1.496 

+1.123 x 10 -3 

1.497 

+5.502 x 10~ 4 

4 

2.628 

-1.380 x 10 -2 

2.629 

-1.285 x 10 -2 

2.640 

-1.260 x 10 -2 

5 

4.023 

-2.758 x 10 -2 

4.037 

-2.147 x 10 -2 

4.075 

-1.942 x 10 -2 

6 

5.554 

-2.543 x 10" 2 

5.562 

-1.935 x 10 -2 

5.632 

-1.610 x 10" 2 

7 

7.144 

-2.618 x 10“ 2 

7.149 

-1.614 x 10“ 2 

7.249 

-1.159 x 10“ 2 

8 

8.759 

-3.077 x 10" 2 

8.767 

-1.745 x 10“ 2 

8.898 

-8.692 x 10" 3 

9 

10.389 

-2.555 x 10“ 2 

10.406 

-1.234 x 10" 2 

10.568 

-7.080 x 10“ 3 

10 

12.033 

-2.250 x 10“ 2 

12.057 

-1.056 x 10“ 2 

12.248 

-6.936 x 10" 3 



(3 = 5 


II 

CO 


(3=1 

mode 

Re(w)/fl 

Im(w)/0 

Re(w)/f2 

Im(w)/Q 

Re(w)/f2 

Im(w)/Q 

1 

0.000 

-4.422 x 10” 1 

0.000 

-6.030 x 10” 1 

0.965 

-2.822 x 10" 1 

2 

0.737 

-1.920 x 10" 1 

0.772 

-1.980 x 10“ 4 

1.501 

-1.595 x 10“ 3 

3 

1.498 

+1.309 x 10 -4 

1.499 

-2.745 x 10 -4 

2.731 

-8.620 x 10 -3 

4 

2.652 

-1.211 x 10 -2 

2.667 

-1.135 x 10" 2 

4.363 

-7.434 x 10 -3 

5 

4.113 

-1.803 x 10 -2 

4.160 

-1.595 x 10" 2 

6.136 

-2.720 x 10 -3 

6 

5.701 

-1.381 x 10 -2 

5.784 

-1.103 x 10" 2 

- 

- 

7 

7.348 

-8.594 x 10 -3 

7.469 

-4.305 x 10" 3 

- 

- 

8 

9.032 

-6.711 x 10" 3 

9.185 

-2.242 x 10" 3 

- 

- 

9 

10.727 

-3.669 x 10" 3 

10.920 

-1.430 x 10" 3 

- 

- 

10 

12.437 

-2.614 x 10" 3 

- 

- 

- 

- 
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Table 3. Eigenfrequencies of One-dimensional Perturbations for T = 0.10 




j3 = oo 


/3 = 100 


o 

T— 1 

II 

mode 

Re(w)/f2 

Im(w)/f2 

Re(w)/f2 

Im(w)/fi 

Re(w)/r2 

Im(w)/fi 

1 

0.354 

-6.363 x 10" 1 

0.363 

-7.086 x 10" 1 

0.379 

-6.983 x 10” 1 

2 

0.809 

-3.897 x 10 -1 

0.817 

-3.963 x 10" 1 

0.855 

-4.454 x 10” 1 

3 

1.614 

-1.078 x 10 -2 

1.617 

-1.173 x 10 -3 

1.619 

-1.517 x 10 -2 

4 

2.740 

+7.693 x 10 -3 

2.741 

+8.745 x 10" 3 

2.759 

+8.427 x 10 -3 

5 

4.316 

-3.565 x 10 -2 

4.319 

-2.545 x 10" 2 

4.377 

-2.297 x 10 -2 

6 

6.019 

-5.747 x 10" 2 

6.052 

-3.498 x 10“ 2 

6.155 

-3.220 x 10" 2 

7 

7.846 

-5.831 x 10 -2 

7.853 

-3.778 x 10“ 2 

7.996 

-2.700 x 10" 2 

8 

9.663 

-8.771 x 10" 2 

9.678 

-3.491 x 10“ 2 

9.866 

-2.868 x 10" 2 

9 

11.509 

-7.421 x 10“ 2 

11.522 

-3.069 x 10" 2 

11.752 

-2.190 x 10“ 2 

10 

13.352 

-8.511 x 10“ 2 

13.375 

-2.640 x 10“ 2 

13.645 

-1.699 x 10" 2 



P = 5 


II 

CO 


P=l 

mode 

Re(w)/fl 

Im(w)/0 

Re(w)/f2 

Im(w)/Q 

Re(w)/f2 

Im(w)/0 

1 

0.894 

-5.006 x 10” 1 

0.935 

-5.721 x 10” 1 

1.033 

-8.362 x 10" 1 

2 

1.622 

-2.052 x 10“ 2 

1.624 

-2.678 x 10“ 2 

1.629 

-5.097 x 10“ 2 

3 

2.776 

+8.532 x 10 -3 

2.796 

+8.896 x 10" 3 

2.877 

+9.385 x 10 -3 

4 

4.433 

-2.016 x 10 -2 

4.498 

-1.702 x 10" 2 

4.745 

-7.109 x 10 -3 

5 

6.250 

-2.726 x 10 -2 

6.361 

-2.199 x 10" 2 

6.776 

-7.802 x 10 -3 

6 

8.132 

-2.718 x 10 -2 

8.286 

-2.087 x 10" 2 

8.863 

-2.309 x 10 -3 

7 

10.042 

-1.719 x 10 -2 

10.241 

-1.584 x 10" 2 

- 

- 

8 

11.964 

-1.949 x 10" 2 

12.207 

-7.872 x 10" 3 

- 

- 

9 

13.896 

-7.560 x 10" 3 

14.180 

-2.576 x 10" 3 

- 

- 

10 

- 

- 

- 

- 

- 

- 
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Fig. 1.— One-dimensional steady-state shock profiles for T = 5% (black) and 10% (red) with 
differing f3. Each dot marks the magnetosonic point. The shock becomes weaker for smaller T and 
/3. Note that shear reversal in the immediate postshock regions is stronger for larger T and (3. 
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Fig. 2.— Five odd-mode eigenfunctions Si for ID perturbations with k y = 0 when (left) (5 = oo 
and (right) f3 = 10. The spiral forcing is set to T = 5%. Red and blue lines represent the real 
and imaginary parts, respectively. All values are normalized such that Re(Si) = Im(Si) = 1 at the 
magnetosonic point located at x/L = 0.007 and x/L = 0.017 for f3 = oo and 10, respectively. The 
vertical line in each panel indicates the shock front. 
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Fig. 3.— (a) Temporal variations of the gas surface density X at x = 0 from a ID simulation with 
J 7 = 10% and /3 = 1. Red solid lines envelope the fluctuation amplitudes of X. The inset zooms 
in the time range 220 < tQ < 260 to display the density fluctuations, (b) The power spectrum of 
the density fluctuations. The frequencies marked by the red arrows represent the real parts of the 
eigenvalues listed in Table 3. 
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Fig. 4.— Non-axisymmetric dispersion relations of the ten lowest-frequency eigenmodes for T = 5% 
and (3 = 10. The modes are numbered in the increasing order of Re(u;) at k y = 0. In each panel, 
the blue solid line (left y-axis) gives Im(cj), while the red dashed line (right y-axis) is for Re(cu). 
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k y L kyL 

Fig. 5.— Dependence on (3 and k y of the growth rates Im(u;) of the dominant overstable modes for 
T = 5%. For j3 > 100, the WI is dominated by the n = 7 mode, while several modes shown have 
similar growth rates for j3 < 10. Star symbols mark the growth rates and wavelengths of the WI 
measured from direct numerical simulations in Section 4.3. 
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Fig. 6.— Dependence on f3 and k y of the growth rates Im(u;) of the dominant overstable modes for 
J- = 10%. For /3 > 5, the WI is dominated by the n = 4 mode, while several different modes have 
similar growth rates for /3 < 3. 
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Fig. 7.— Eigenfunctions Si, U\, V\, M\, and Ej of an overstable mode (left) with cn/fl = 84.4+1.13* 
and k y L = 78.1 for /3 = 100 and T = 5% and (right) with lj/Q = 59.1 + 0.50* and k y L = 53.4 for 
/3 = 10 and T = 5%. The red solid and blue dotted curves represent the absolute values of the real 
and imaginary parts. The shock front is indicated by the vertical dashed line in each panel. The 
black dots in the top panels mark the magnetosonic points. The black solid lines enveloping the 
eigenfunctions in the bottom panels draw the solutions of Equation (A9). 
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Fig. 8.— Distributions of the perturbed magnetic fields (black solid lines) in the regions with 
—0.2 < x/L < 0.3 and 0 < yj\ y < 1 overlaid over the perturbed PV, Re(£i), displayed in color 
scale for the (a) j3 = 100 and (b) [3 = 10 cases with T = 5%. White dots in both panels trace the 
wavefronts of the perturbed PV. Colorbar labels Re(^i)/(DS“ 1 ). 
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Fig. 9.— Ratio of the magnetic to hydrodynamic terms in the expression for the PV jumps across 
a disturbed shock front due to (a) the perpendicular velocity U\ and (b) the distortion amplitude 
Z\. Filled and open circles correspond to T = 5% and 10%, respectively. 
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Fig. 10. — Snapshots of the gas surface density as well as the configuration of the magnetic fields 
in the regions with —0.2 < x/L < 0.35 and 0.5 < y/L < 1.2 from 2D simulations with (3 = 100, 10, 
5, and 3 at if) = 8, 16, 25, and 40 from left to right, respectively. The arm forcing is T = 5% and 
the grid resolution over the domain size of L x 2L is 2048 x 4098. The number of the nonlinear 
features grown by the WI along the y-direction is 21, 15, 12, and 9 from left to right. Colorbar 
labels log(E/£ c ). 
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Fig. 11.— Temporal evolution of the maximum surface density measured at x = 0 from the models 
shown in Figure 10. The growth rate measured from the slope indicated as the line segment in each 
model is consistent with the results of the normal-mode linear stability analysis. 
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Fig. 12.— Difference between the pitch angles of the main arm and the nonlinear structures 
stretched from it in the numerical simulations with T = 5%, averaged over the regions with 
x s h < x < x s h + L/2. Filled circles and errorbars give the mean values and standard deviations 
over the time interval of At = 5/D from the time epoch shown in Figure 10. 



